# This code checks the sparsity pattern of the inverse spectral density matrix 
# for a multivariate autoregressive sequence. It uses the notation and conditions
# in Songsiri et al. There are other ways to check but this has the benefit that it
# produces the matrices required in that paper. 
# 
# inputs: 
#   A = [A_1,...,A_p], autoregressive coefficients.
#   C, the noise covariance matrix
#   p, the problem dimension. 
#   k, the lag 
#   eps, a threshold for sparsity.
#
# returns:
#   M, the p(k+1) by p(k+1) matrix [I,A]^T C^(-1) [I,A] as in Songsiri et al. 
#   D,  the p by p(k+1) matrix D(M) with notation as in Songsiri et al. 
#   Sparsity, a p by p matrix encoding the sparsity. See Songsiri et al. for conditions. 
##########################################################################################

check_sparsity <- function(A, C, p, k, eps = 10^-6){

  # construct M
  M <- t(cbind(diag(p), A)) %*% solve(C, cbind(diag(p), A))
   
  # construct D and recover Sparsity
  # construct D_0
  D <- matrix(0, p, p)
  Sparsity <- matrix(0, p, p)
  for( index.block in 0:k ){
    D <- D + M[(index.block * p + 1):((index.block + 1) * p), (index.block * p + 1):((index.block + 1) * p)]
  }
  
  Sparsity <- (abs(D) > eps)
  
  # consruct D_i (which is a sum of matrix block off-diagionals) for i in 1:k
  for(i in 1:k){
    D.temp <- matrix(0,p,p)
    for( index.block in 0:(k-i) ){
      block.rows    <- ((i+index.block)*p+1):((i+index.block+1)*p) 
      block.columns    <- (index.block * p + 1):((index.block + 1) * p)
      D.temp <- D.temp + 2 * M[block.columns, block.rows ]
    }
    D <- cbind(D, D.temp)
    Sparsity <- Sparsity & (abs(D.temp) > eps)
  }
  
  storage.mode(Sparsity) <- "integer"
  list("M" = M, "D" = D, "Sparsity" = Sparsity)
}
